/* 
   Example of a hybrid API using OpenMP and MPI 
   . Sums the elements of a 1D array of size n.
   . Uses openMP reduction and MPI_Reduce()
	 
   using a gcc-based MPI
   > mpic++ -fopenmp reduce.cc -o reduc
   using a icc-based MPI
   > mpic++ -openmp reduce.cc -o reduc
   and run with
   > qsub -I -l nodes=2
   > export OMP_NUM_THREADS=4; mpirun -np 2 ./reduce [array_size]
   The idea here is to split the entire array into MPI_NUM_PROCS * OMP_NUM_THREADS subintervals and then pick the interval to be the interval passed to the next iteration
*/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <iostream>
#include <omp.h>
#include <mpi.h>
#include <vector>
#include <unistd.h>

using namespace std;

#define N_DEFAULT 20
#define REAL float
#define USE_GPU 0
#define MAIN_RETURN_NOERROR 0

/*
Imitate the interface of bsearch
void *bsearch(const void *key, const void *base,
	      size_t nmemb, size_t size,
	      int (*compar)(const void *, const void *));
*/

int userComparOperator(const void *input1, const void *input2);

int getSplitters(int* splitter_bnds, int *splitters, int size_nmemb, int mpi_size);

void *pbsearch(const void *key, 
	       const void *base, 
	       size_t size_nmemb, 
	       size_t size,
	       int (*compar)(const void *, const void *),
	       MPI_Comm comm,
	       int ntds);

//#define ucmpop(x1,x2) userComparOperator( (const void *) (x1), (const void *) (x2) );

	       
int main( int argc, char *argv[]){

  int n; // original problem size
  n= (argc==2) ? atoi(argv[1]) : N_DEFAULT;	
  int m = n; // interval size at each iteration
  int namelen;

  
  int i;
  int r = 0;

  double key = 0.3000;
  int ntds = 2;	  // number of threads for each process

  MPI_Init(&argc,&argv); 	// intialize MPI and MPI-related variables 
  MPI_Comm comm=MPI_COMM_WORLD;

  vector<double> A;  A.resize(n);

  // initialize array values
  for(int i=0;i<n;i++) A[i] =(1.0 / n) * i;  
      //	    cout << A[i] << endl;

  //double *bucket_begin = A.begin();
  int bucket_size = n;
  //  printf("hello\n");
  pbsearch(&key, 
	   &A[0],
	   bucket_size, 
	   sizeof(A[0]),
	   userComparOperator,
	   comm,
	   ntds);

  MPI_Finalize();

  return 0;
}

int userComparOperator( const void *x1, const void *x2){
	double *a = (double *) x1;
	double *b = (double *) x2;

	if (*a > *b) 	  return 1;
	else if(*a < *b) return -1;
	else return 0;
}

void *pbsearch
(
 const void *key,      // address of the key on each process
 const void *base,     // address of the base on each process
 size_t size_nmemb,    // number of the elements of the array
 size_t size,          // size of array element
 int (*compar)(const void *, const void *), //compare funciton pointer
 MPI_Comm comm,        // MPI communicator
 int ntds              // Number of OpenMP threads 
)              
{
  
  char name[32];
  int mpi_rank; MPI_Comm_rank( comm, &mpi_rank);
  int mpi_size; MPI_Comm_size( comm, &mpi_size);

  //  int bucket_size = int(ceil(double(m)/ntds));
  //  int nblks =  int(ceil( double(m)/blk_size));

  //MPI_Get_processor_name(name,  &namelen);
  int root = 0;
  bool iamroot = (mpi_rank==root);
  int *splitter_vec;
  int splitters[2*mpi_size];

  // Compute the splitters and store them in a array of integers
  if (iamroot) 
    {

      int splitter_bnds[2];
      splitter_bnds[0] = 0;
      splitter_bnds[1] = size_nmemb;
      getSplitters(splitter_bnds, splitters, size_nmemb, mpi_size);
      splitter_vec = splitters;  
      for (int j = 0; j < 2*mpi_size; j++)
	cout << splitters[j] << endl;
    }

  int mysplitter[2];
 
  //  MPI_Gather(results, 2, MPI_INT, myresult, 2, MPI_INT, 0, comm);
  omp_bsearch_seg(base, key, splitters, num_elem, size_elem, op, ntds); 
  printf("my rank % i: %i %i\n", mpi_rank, mysplitter[0], mysplitter[1] );

  // Here we can wrap up the following codes to be a function
  /* 
     Input: (
      - address of array A
      - splitter of array A
      - the comparison operator
      - number of elements

     Output:
      - flag of whether of A is found
      - if found pass back the index of the array, otherwise passback -1

     prototype:
      - 
      omp_bsearch
      
  */

  cout << "iterations " << iter << endl << endl;
  cout << "left pivot" << *glsplt  << "      "
       << "right pivot" << *grsplt << endl;

  return 0;
}


int getSplitters(int* splitter_bnds, int *splitters, int size_nmemb, int mpi_size)
{

  int bucket_size = int( floor( double(size_nmemb)/mpi_size));
  int bigger_size = size_nmemb%mpi_size;
  int left_splitter = splitter_bnds[0];
  int right_splitter = splitter_bnds[0];

  for (int i = 0; i < bigger_size; i++)
    {
      left_splitter = right_splitter;
      right_splitter = left_splitter + bucket_size + 1;
      splitters[2*i] = left_splitter;
      splitters[2*i+1] = right_splitter;
    }

  for (int i = bigger_size ; i < mpi_size; i++)
    {
      left_splitter = right_splitter;
      right_splitter = left_splitter + bucket_size;
      splitters[2*i] = left_splitter;
      splitters[2*i+1] = right_splitter;
    }
  
  return 0;
  


}
 


//  omp_bsearch_seg(base, key, splitters, num_elem, size_elem, op, ntds); 	
int omp_bsearch_seg
(
 void* base,
 void* key,
 int num_elem,
 int size_elem,
 int (*compar)(const void *, const void *), //compare funciton pointer
 int ntds                                   // Number of OpenMP threads 
)
{	
#pragma omp parallel
#pragma omp master
  omp_set_num_threads( ntds );

  //to make it generic we need to make them void instead
  double *glsplt = (double *)base + mysplitter[0];
  //double *glsplt = (double *)base; 

  double *glsplt_new;
  double *lsplt;
  double *grsplt = (double *)base + mysplitter[1];
  //double *grsplt = (double *)base + size_nmemb;
  double *grsplt_new;
  double *rsplt;

  int iter = 0;
  int blk_size=1000;
  int nblks;
  int lid, rid;
  int m = 10000;
  int tdn;

  while( blk_size > 2 )
    {

      m = grsplt - glsplt + 1;

      blk_size = int(ceil(double(m)/ntds));
      nblks =  int(ceil( double(m)/blk_size));

      /*
      cout << "m"  << m << endl;
      cout << "nblks" << nblks << endl;
      cout << "blk_size" << blk_size <<endl;
      */

      // inner loop of parallel for loop
#pragma omp parallel for private(lsplt, rsplt, tdn, lid, rid) \
                         shared(glsplt, grsplt, glsplt_new, grsplt_new)
      for(int i = 0; i < nblks; i++)
	{
	  tdn = omp_get_thread_num();
	  printf("%f , %f\n", *glsplt, *grsplt);

	  lsplt = glsplt + i * blk_size;

	  if (i == nblks - 1)
	    rsplt = grsplt;
	  else
	    rsplt = glsplt + (i+1) * blk_size;

	  // We only need to pass the pointer of the key and the left/right splitter 
	  lid = compar(lsplt, key);
	  // lid = ( *lsplt >= key);
	  rid = compar(rsplt, key);
	  // rid = ( *rsplt > key);

	  printf("from thread %i, left spliter %f (%i), right spliter %f(%i)\n", tdn, *lsplt, lid, *rsplt, rid);
	
	  if ( (lid == 0 || lid == -1) && rid == 1)
	    {
	      glsplt_new = lsplt;
	      grsplt_new = rsplt;
	    }

	}

      // Update global splitters
      glsplt = glsplt_new;
      grsplt = grsplt_new;

      // cout << *grsplt << " " << *glsplt << endl;
      iter ++;

    }
}
